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Abstract 

Given a Feynman parameter integral, depending on a single discrete variable N and a real 
parameter e, we discuss a new algorithmic framework to compute the first coefficients of its 
Laurent series expansion in e. In a first step, the integrals are expressed by hypergeometric multi- 
sums by means of symbolic transformations. Given this sum format, we develop new summation 
tools to extract the first coefficients of its series expansion whenever they are expressible in terms 
of indefinite nested product-sum expressions. In particular, we enhance the known multi-sum 
algorithms to derive recurrences for sums with complicated boundary conditions, and we present 
new algorithms to find formal Laurent series solutions of a given recurrence relation. 
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Introduction 



Sta r ting with s i ngle summation oyer hypergeometric terms developed, e.g . . mltjosperl 
(1978); Zcilbcrgcr ( ll990al) ; |Petkovsekl(|l992l) ; lAbramov and Petkovsel1(|l994l ): |Paulel(|l995l) 



Gosperl 



symbolic summation has been intensively enhanc e d to multi-summation l i ke, e.g., the 
holono mic approach of IZeilbergerl (Tl990bh : IChv zakl (l200dh ; ISchneiderl (l2005al) ; iKoutschanl 
20091) . In this article we use the techniques of iFasenmverl (Il945l); IWilf and Zeilberge 



19921) which lead to efficient algorithms developed, e.g., in IWegschaiderl (|l997t ) to com 



pute recurrence relations for hyper geometric multi- sums. Besides this, we rely on multi- 
summation algorithm s presented i n Schneider ( 20071 ) that generalize the summation tech- 
niques worked out in IPetkovsek et al. ( 19961); the underlying algo rithms are based on a 
refined difference field theory elaborated in ISchneiderl (|2008l l2010h that is adapted from 
Karr's IlS-ficlds originally introduced in iKarrl (Il98lh . 

We aim at combining these summation approaches which leads to a new framework 
to evaluate Feynman integrals. In a nutshell, given a Feynman integral, we transform 
it to hypergeometric multisums, compute afterwards linear recurrences for these multi- 
sums, and finally decide constructively by recurrence solving whether the integrals (resp. 
the multisums) have series expansions whose coefficients can be represented in terms 
of indefinite nested sums and products. The method consists of a completely algebraic 
algorithm. It is therefore well-suited for implementation in computer algebra systems. 

We show in a first step that Feynman parameter integrals, which contain local operator 
insertions, in Z?-dimensional Minkowski space with one time- and (D — 1) Euclidean space 
dimensions, e = D — 4 and e 6 R with |e| -C 1, can be transformed by means of symbolic 
computation to hypergeometric multi-sums S(e, N) with N an integer parameter. Given 
this representation, one can check by analytic arguments whether the integrals can be 
expanded in a Laurent series w.r.t. the parameter e, and we seek summation algorithms 
to compute the first few coefficients of this expansion whenever they are representable 
in terms of indefinite nested sums and products. Due to the special input class of Feyn- 
ma n integrals, these solutions c an be usually trans f ormed to ha r moni c sums or S- sums ; 
Blumlein and Kurthl (|l999i ); IVermaserenl (|l999f ); lMoch et all (|2002l) ; lAblingerl (|2009f ). 



see 



In general, we present an algorithm (see Theorem 1) that decides constructively, if 
these first coefficients of the e-expansion can be written in such indefinite nested product- 
sum expressions. Here one first computes a homogeneous recurrence by WZ-thcory and 
Wegschaider's approach. This recurrence together with initial values gives an alternative 
representation for the series expansion (see Lemma 1). Moreover, we develop a recur- 
rence solver (see Corollary 1) which computes the coefficients of the expansion in terms 
of indefinite nested product-sum exp ressions wheneve r this is possible. The backbone o f 
this solver relies on alg orithms from Petkovsek ( 19921 ) ; lAbramov and Petkovsek ( 1994 ) ; 
Schneider! (|200ll , l2005bl ). Since the solutions are highly nested by const ruction, their sim - 
plification to sum representations with minimal depth are crucial; see Schneider ( 20101 ). 

From the practical point of view there is one crucial drawback of the proposed so- 
lution: looking for such recurrences is extremely expensive. For our examples arising 
from particle physics the proposed algorithm is not applicable considering the available 
computer and time resources. On that score we relax this very restrictive requirement 
and search for possibly inhomogeneous recurrence relations. However, the input sums 
have summands which present poles outside the given summation ranges . Combining 
Wegschaider's package MultiSum and the new package FSums presented in IStanl ( 20101 ) 
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we determine recurrences with inhomogeneous sides consisting of well-defined sums with 
fewer sum quantifiers. Applying our method to these simpler sums by recursion will lead 
to an expansion of the right hand side of the starting recurrence. Finally, we compute 
the coefficients of the original input sum by our new recurrence solver mentioned above. 

The outline of the article is as follows. In Section 2 we explain all computation steps 
that lead from Feynman integrals to hypcrgeomctric multi-sums of the form (7) which can 
be expanded in a Laurent expansion (If) where the coefficients Fi(N) can be represented 
in the form (12). In the beginning of Section 3 we face the problem that the multi- 
sums (7) have to be split further in the form (13) to fit the input class of our summation 
algorithms. We first discuss convergent sums only. The treatment of those sums which 
diverge in this special format or sums with several infinite summations that have difficult 
convergence properties will be dealt with later, cf. Remark 5. In the remaining parts of 
Section 3 we present the general mechanisms to compute the first coefficients Fi(N) for 
a given hypergeometric multi-sum. In Section 4 we present an algorithmic approach to 
hypcrgeometric sums with non-standard boundary conditions. This allows us to generate 
the inhomogeneous sides of recurrences delivered by Wegschaider's package MultiSum. 
Finally, in Section 5 we obtain a method that is capable of computing the coefficients 
Fi(N) in reasonable time. Conclusions are given in Section 6. 



2. Multiple sum representations of Feynman integrals 

We show how integrals emerging in renormalizablc Quantum Field Theories, like Quan- 
tum Electrodynamics or Quantum Chromodynamics, see e.g. iBliimlein (l2009h . can be 



transformed by means of symbolic computation to hypcrgeomctric multi-sums. We study 
a very general class of Feynman integrals which are of relevance for many physical pro- 
cesses at high energy colliders, such as the Large Hadron Collider, LHC, and others. 
The processes obey special-relativistic kinematics with energy-momentum vectors in 



Minkowski space, M , see e.g.. iNaas and Schmidl ([19611 ). i.e., a D-dimcnsional linear 



space where the elements a = (do, a) £ M. D decompose into the time coordinate a,Q £ 
K and the spatial coordinates a £ M. ^ 1 which form a D — 1-dimcnsional Euclidean 
subspace; the bilinear form is defined by a. b = (a, b) = aobo — ab £ K for b = (bo, b) £ M D . 
Below analytic continuations in D := 4 + e with e £ K are considered. Here we study 
integrals 

1 ,PJ J (2ir) D J (2tt) d {-p\ + m\) l i . . . {-p 2 + m|)'* V w 

with A,p,pi £ M D and m l £ {0, M} for some M £ R with M > 0. The restriction that 
there is only one mass M is the only one specifying the class of Feynman diagrams from 
arbitrary ones. The propagator powers U obey U £ N and for the special vector A in 
(1) one has A. A = 0. The numerator AT is usually given in terms of finite sums where 
the range depends on a discrete parameter N and where the summand depends on the 
scalars p-Pj,Pi-Pj, A.pi (1 < i, j < k), on M 2 and on N. In particular, for each N £ N, 
Af is a polynomial in terms of these scalars and M 2 where the exponents of the A.pj 
(1 < j < k) in a given monomial sum up to N and the exponents of the remaining scalars 
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and M 2 are constant. The Sy occurring in (1) are shortcuts for Dirac delta functions in 
D dimensions 5v = (Ym=i a v,Wi) , av,i S Q. I.e., if ay,i ^ 0, we get 

k 

L 




d D Pt 5^ [y a v , lPl )f( Pi ):={^\ with u := -— V ov.jPU (2) 



here / stands for the integrand of (1). For each such rule (2) for the remaining Sy, one 
integral sign in (1) can be eliminated. As a consequence we obtain integrals of the same 
shape but with fewer integral signs. Suc h an integral ma y be easily lin early transformed 
into Euclidean integrals (Wick rotation, Fevnman ( 19491 ); Iwick ( 1950( )) in the Euclidean 
space by replacing a = (a,Q,d) £ M D with a = (iao,a). In this way, for b = (bo,b) the 
bilinear form (a, b) — — ao&o — a.b < obtains a definite sign; \J — {a, a) is then the 
Euclidean norm ||a||. Summarizing, we obtain an Euclidean integral of the same shape 
as (1) with the Euclidean momenta pi,p (instead of Pi,p) and where the denominators 
can be written in the form ((J2j=i c ^Pj) 2 + m 2 ) li with cy g Q (instead of (— p 2 + m 2 ) li ); 
this format is due to the usage of (2). 

Subsequently, we show how this Euclidean integral can be mapped to an integral on 
an m-dimensional unit cube. Define D{ := (52j=i c ^Pj) 2 + TO ?- Then we loop over r 
(r = 1,2,..., A) as follows. For the rth iteration, fix q := p r . W.l.o.g. assume that 
c£ G {0, 1} for 1 < i < k. Now collect those denominator factors where q occurs, 
say rij=i ^i- ( n e f^)- Then we use the formula 

if? - rpw I ' ^ • • • L ' dx -i ti 11 - w A^.r (3) 

with / = YTj—i hj ; here 5 is the Dirac delta function, the variables Xk are called Feynman 
parameters, and T(z) denotes the Gamma-function. Due to the Dirac delta function, we 
get that A := XiD^ + . . . x n D ij = q 2 + a.q + b where a and b are expressions free of q. 
Hence we can write A = (q + a/2) 2 + R with R := —a 2 /4 + b being free of q. Replacing the 
denominator of our integral by this formula, we can simplify A further. Namely, using 
the shift-invariance w.r.t. the vector q, which holds in D-dimensional Euclidean space, 
the denominator A can be brought to the form (g 2 + R) without changing the integral. 
Finally, expanding the numerators and applying the g-integral termwise lead to integrals 

of the form J ^^qr^p ~ where the expression q\ is free of q. If m is odd, i.e., an odd 
number of vector multiplications w.r.t. q arise, the integral evaluates to by symmetry. 
If m is even, one exploits the simplification 

d°q Yl™l\qx-q _ r{D) f d D k (fy 



{2tt) d (q 2 +R) 1 K 'J (2ir) D (q 2 +RY 
where r(D) stands for a rational function in D (i.e., in e) that can be determined by an 
explicit formula. To this end, the following formula is applied to the remaining integrals: 

d D q (fY _ 1 r(r + D/2)T(l -r - D/2) 
(2tt) d {q 2 + R) 1 ~ (16tt 2 ) b / 4 T(D/2)T(l)(R 2 y- r - D / 2 ' 

Usually, these operations are carried out in terms of tensors to keep the size compact 
and to determine additional relations efficiently. The above procedure is repeated until 
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all momentum integrals for the p r (r = 1, 2, . . . , k) are computed. As a result one is left 
with the integrals over x, G [0, 1], equipped with a pre-factor C(e,N,M). 

Step 1: From Feynman parameter integrals to Mellin-Barnes integrals and multinomial series. 

Parts of these scalar integrals again can be computed trivially related to the ^-distributions, 



dxi 8 



fe=i 



x k - 1 



1 



e x k ) n ^™)> 

k=l,kjtl ' m=l,m#i 

where 0(.z) is 1 if z > and otherwise. There may be more integrals, which can be 
computed, usually as indefinite integrals, without special effort. Mapping all Feynman- 
paramctcr integrals onto the m-dimcnsional unit cube (as described above) one obtains 
the following structure : 



X(e,N) = C(e,N,M) / dy 



dy, 



i,i(e,N) 



(4) 



with k G N, r\ , . . . , r k G N and where 
/3(e) G 0(e), and similarly onAe, N) = n 



(e) is given by a rational function in e, i.e., 
.,.N + cii,i for some n it i G {0, 1} and a^i G Q(e), 
see also iBogner and Weinzierl ( 2010T ) in the case when no local operator insertions are 
present. C(e, N, M) is a factor which depends on the dimensional parameter e, the integer 
parameter iV and M. Pi(y), Q(y) arc polynomials in the remaining Feynman parameters 
y = (j/i, . . . , y m ) written in multi-index notation. In ( 4) all terms which s tem fro m local 
operator insertions were geometrically resummed; see iBierenbaum et al. | (l2009bl) . 

Remark. (1) After splitting the integral (4) (in particular, the k summands), the inte- 
grands fit into the input class of the multivariate Almkvist-Zeilberger algorithm. Hence, 
if the split integrals are properly def ined, they obey homogen e ous re currence relations in 
N due to the existence theorems in Apagodu and Zeilberger (|2006h . However, so far we 
failed to compute these recurrences due to time and space limitations. 
Remark. (2) Usually the calculation of X(e, N) for fixed integer values of AT is a simpler 
task. If sufficiently many of these values arc known, one may guess these recurrences 
and with this i nput derive closed forms for X(e,N) using the techniques applied in 
Bliimlein et al. I (|2009h . This has been illustrated for a large class of 3-loop quantities. 
However, at present no method is known to calculate the amount of moments needed. 

The y^-integrals finally turn into Euler integrals. Here we outline a gen eral framework 



although in pr actice, different algorithms arc used in specific cases, cf. e.g. lAblinger et al 



( 2010a . 2011bl ). To compute the integrals (4) over the variables j/j we pr oceed as follows: 

• decom pose the denominator function using Mellin-Barnes integrals, see lParis and Kaminski 
(|2001l ) and references therein, 

• decompose the numerator functions, if needed, into multinomial series. 
The denominator function has the structure 



lQ(y)} 



^2ik(y) 



Lk=l 



with q k {y) = Oi . . - dm where G 1 — J/i} for 1 < i < m. This function can be 

decomposed applying its Mcllin-Barncs integral representation (n — 1) times, 



(A + B)i 



1 

2ni 



da A a B~ q - a 



T(-a)T(g + a) 



(5) 
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Here 7 denotes the real part of the co ntour. Often Eg. (5) h a s to be considered in 
the sense of its analytic continuation, see Whittaker and Watson ( 19961 ). The numerator 
factors [Pi,i(y)] ai ' ,< - S ' N) obey 



[PiAv)] 



.k=l 



a*,l(e,JV) 



where the monomials p k (y) have the same properties as q k (y) ■ One expands 



[Pi,i(y)] ai ' ,(e ' N) = £ 



a h i{e,N) 
fci, . . . , k w —\ 



w — l 



\{vi{y) kl P w {yr Ae ' N) ~^' 



1=1 



Now all integrals over the variables yj can be computed by using the formula 



<W _1 (l-») 



p-i 



B(a,/3) 



r(q)r(/3) 
r(a + P) 



and one obtains 



I(e, N) 



(2m) n 



71+ioo 



£ 

fcl=l 



£„(W,ki,...,A„-i) J 



£ £c fe ( £ ,iv,M) 



fc=i 



r(zi ;fc ) . . .r(z M|fc ) 



(6) 



£ G N and the summation over fcj comes from the multinomial sums, i.e., the upper 
bounds Li(N), . . . , L V (N, ki, . . . , k v -i) are integer linear in the dependent parameters or 
00. Moreover, the z u & are linear functions with rational coefficients in terms of e, the 
Mcllin-Barncs integration variables <7i , . . . , cr n , and the summation variables ki , . . . , k v . 

Step 2: Representation in multi-sums. The Mcllin-Barncs integrals are carried out applying 
the residue theorem in Eq. (6). The following representation is obtained: 



00 00 L\(N) 

iMD=£~£ E ■ 

n 1 — 1 n r — 1 fci=l 



, £ J2c k (e,N,M) -^hL_ r M v ( 7 ) 



fc»=i 



fe=i 



r(i-u'+i,fc) ■ • ■ r(^io',/c) 



Here the t;.fe are linear functions with rational coefficients in terms of the m,... ,n r , of 
the k\, . . . , k v , and of e. Note that the residue theorem may imply more than one infinite 
sum per Mellin-Barnes integral, i.e., r > n. 

In general, this approach leads to a highly nested multi-sum. Fixing the loop order of the 
Fcynman integrals and restricting to certain special situations usually enables one to find 
sum re presentations with fewer summation signs. E.g., as worked out in lBierenbaum et al 



(|2008f ). one can identify the underlying sums in terms of generalized hypergeometric 
functions, i.e., the number of infinite sums are reduced to one or in some cases to zero. 



Step 3: Laurent series in e. Eq. (7) can now be expanded in the parameter e using 

r(n)r(l + e) 



with £ 



for some r G 



T(n + 1 + s) 
and 

(-£) 



B(n, 1 + e) 



B{n, 1 +£) = - cxp Mr K -^S k (n) = - £(-£")%, . . . ,i(n) 



(8) 
(9) 



\k=l 



k=0 
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and other well-known trans fo rmations for theT- functions. Here the harmonic sums Sg(N) 
Blumlein and Kurth (|l999t) : IVermaserenl f|l999h for TV £ N are recursively defined by 



S b ,s(N) 



N 

E 

k=l 



(sign(6)) fc 



(10) 



Note that in (8) n may stand for a linear combination of parameters with coefficients 
in Q. In case of non-integer weight factors r, f or the par a meters in n analytic continua - 
tions of harmonic sums have to be considered Bliimlcin (|2000L l2009h : Bliimlcin (|2010ft : 
Blumlein and Mochl (|2005l) . In case that n is not an integer one may shift to n — > k n G N, 



which leads to the usual definition of the harmonic sums in (9). However, the summation 
operators have now to be generalized and o ne usually ends up with cyclotomic harmonic 
sums worked out in lAblinger etaH (|2011al ). 

Applying (8) with (9) to each factor in (7) produces for some L > the expansion 



l{e,N) 



OO 

E 

1—-L 



(11) 



L equals the loop order in case of infra-red finite integrals; otherwise, L may be larger. 

Remark 1. In order to guarantee correctness of this construction, i.e., performing the 
expansion first on the summand level of (7) and afterwards applying the summation on 
the coefficients of the summand expansion (i.e., exchanging the differential operator D s 
and the summation quantifiers) analytic arguments have to be considered. For all our 
computations this construction was possible. 

The general expression of the functions Ii (N) in terms of nested sums are 



oo il(iV) 

E E 



,fc„-i 



^2 Y^ H j( N ;ni,...,n r ;,ki,...,k v ) 



(12) 



rii, n r] i ki, k v )y, 

Hj{N; ni, •••> k v ) denote proper hypergeometric termfQ and Sg i . (Lij{N; ni,...,k v )) are 
harmonic sums with the index set o» j and Li t j (usually integer linear) functions of the 
arguments (N; ni,...,k v ). The su m-structure i n (12) is usually obtained performing the 
synchronization of arguments , see Vermaseren ( 19991 ). and applying the associated quasi- 
shuffle algebra, see Blumlein ( 2004) . 



3. First approach to the problem 



In the following we limit the investigation to a sub-class of integrals of the type (1) 
and consider two- and simpler three-loop diagrams, which occur red in the calcu l ation o f 
the massive Wi l son coefficients for deep-inelastic scattering; see Ablinger et al. ( 2011bl ); 



Blumlein et"aTl (|2006h : iBierenbaum et all (|2007t l2009al 120081 ). Looking at the reduction 



1 For a precise definition of proper hypergeometric terms we refer, e.g., to IWegschaideil 1119971 ). For all 
our applications it suffices to know that Hj might be a product of Gamma-functions (occurring in the 
numerator and denominator) with linear dependence on the variables JV, ni, ki times a rational function 
in these variables where the denominator factors linearly. 
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steps of the previous section we obtain the following result. If we succeed in finding 
the representation (11) with (12) it follows constructively that for each N £ N with 
N > A for some A G N the integral I(e,N) has a Laurent expansion in e and thus it is 
an analytic function in e throughout an annu l ar region centered at where the pole at 
e = has order L. In lBierenbaum et all (|2008l ); lAblinger et all (l2011bl . l2010bl ) we started 



with the sum representation of the coefficients (12) and the main task was to simplify 
the expressions in terms of harmonic sums. In this article, we follow a new approach that 
directly attacks the sum representation (7) and searches for the first coefficients of its 
e-expansion (11). By splitting (7) accordingly (and pulling Gk{e, N , m)) our integral can 
be written as a linear combination of hypergeometric multi-sums of the following form. 

Assumption 1. 

oo oo Li(N) L 2 (N,j ) L r (N,jo,--;jr-i) 

S(s,N)= £ ■■■ E E E ••• E F{N,v,j ,...,jr-i),e) (13) 

cri=pi o- 3 =p 3 jo=qo ji=qi j r =q r 

where 

(1) N E N with N > A for some given A E N, e > is a real parameter; 

(2) the upper summation bounds Li(N, jo, ■ • ■ , ji-i) are integer linear in iV, j , . . . ,ji-i, 
and the lower bounds are given constants p il qi E N for all 1 < i < s and < I < r; 

(3) J 7 is a proper hypergeometric term (see Footnote 1) with respect to the integer 
variable N and all summation variables (cr,j) = (o\, . . . ,a s ,jo, ■ ■ ■ ,j r ) € Z s+r+1 . 

Remark 2. While splitting the sum (7) into sums of the form (13) it might happen 
that the infinite sums over individual monomials diverge for fixed values of e, despite the 
convergence of the complete expression. We will deal with these cases in Section 5 and 
consider only sums which are convergent at the moment. 

In other words, we assume that (13) itself is analytic in e throughout an annular region 
centered at and we try to find the first coefficients F t (N), F t +i(N), . . . , F U (N) in terms 
of indefinite nested product-sum expressions of its expansion 

5(e, N) = F t (N)e* + F t+1 (N)e t+1 + F t+2 (N)e t+2 + .... (14) 

with t E Z. In all our computations it turns out that the summand T (N, a, j, e) satisfies 
besides properties (l)-(3) the following asymptotic behavior: 

(4) for all 1 < i < s we have 

J" (JV.tr, j,e) = O (V dle ~ C,<T *) as oW°o with c, > 0, d t > 0; (15) 

for simplicity we do not consider the log-parts. For later considerations in Section 4 
we suppose that such constan ts Cj and di are given explicitly. E.g., using the behavior 
( Whittaker and Watsonl . [l99l Section 13.6) of logT( z) for large \z\ in the region where 



arg(z)| < 7r and |arg(z + a)\ < ir: 

log T(z + a) = (z + a - i) log z - z + 0(1), (16) 

such constants can be easily computed. If not all a > for 1 < i < s, things get more 
complicated and -for simplicity- we restrict ourselves to the case that s = 1 and ci = 0; 
we refer again to Section 5 for further details how one can treat the more general case. 
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(5) If s = 1 and c± = 0, we suppose that we are given a constant oGN such that 

oo 

S{e,N)= ^(N,<j 1:J: e) (17) 

CT1=P1 

converges absolutely for any small nonzero e around 0, N > B and any j that runs 
over the summation range. 
Using, e.g., facts about hypergeometric functions from ( Andrews et al. . 19991 Thm. 2.1.1) 



the maximal such constant o in (17) can be determined. 

Example 1. The following sum is a typical entry from the list of sum representations 
for a class of Fcynman parameter integrals we computed: 

00 jy-3 iv-io-3 io+i /, + ^ ("-ff-^rfr + 32 + 2)r( J1 + n + 3) 



oo JV-3 N-j -3j + l / . 

"MO-EE e e ,o 3 



h ) N\ 

(18) 

(I + 1 U (-£)<TiO'i +J2 +3) gl (3- f)^ r(N-j - 1)T(N- ji - j 2 - 1) 

Cii + 4) CT1 (-§ +ji +32 + 4) cti (4 - |) . i+ . 2 r( ffl + i)r(ii + 4)r(jv - j - 2) ; 

we denote by (x)k = x(x + 1) . . . (x + k — 1) the Pochhammer symbol defined for non- 
negative integers k. Then using formulas such as (x)k = T(x + k)/T{x) and (^) = 
T{x + l)/T{x — k + l)/T{k+l) and applying (16) we get the asymptotic behavior C(trf 5 ) 
of the summand. Moreover, we choose the maximal o = 3 such that condition (17) is 
satisfied. 

Subsequently, we will develop an algorithm that finds, whenever possible, representations 
for the coefficients in the expansion (14) in terms of indefinite nested sums and products- 
Theorem 1. Let S(e,N) be a sum with properties (l)-(5) from Assumption 1 which 
forms an analytic function in e throughout an annular region centered at with the 
Laurent expansion (14) for some t G Z for each nonncgativc N; let ueN. Then there is 
an algorithm which finds the maximal r £ {t — 1, t, . . . ,u} such that the ft(N), . . . , f r (N) 
are expressible in terms of indefinite nested product-sums; it outputs such expressions 
F t (N) , . . . , F r (N) and A £ N s.t. f t (k) = F t {k) for all < i < r and all k £ N with k > A. 

This result is based on the fact that such sums S(e,N) satisfy a recurrence relation. 
Example 2. Consider the single nested sum 

gfe *> - e l ~V k + f ( " " c)r(t + ," )r(JV)r( r i+k+2} P») 

S L(2-|)r(-£ + fc + 4)r(§ + A: + 3)r(A^-fc) V ; 

over a proper hypergeometric term; note that an expansion (14) with t = exists follow- 
ing the arguments from Remark 1. In the first step we compute the recurrence relation 

a (e, N)S(e, N) + oi(e, N)S(e, N + 1) + a 2 {e, N)S{e, N + 2) = h(e, N) (20) 



2 This means in particular indefinite nested sums over hypergeometric terms (like binomials, factorials, 
Pochhammer symbols) that may occur as polynomial expressions with the additional constraint that 
the summation index ij of a sum f(ij) may occur only as the upper index of its inner sums and 

produ cts, but not inside the inner sums themselves; for a formal but lengthy definition sec Schneider 
| |201C]| ), Typical examples are sums of the form (10) above, or of the forms (33) and (34) given below. 
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with 

h(e, N) = -24N - 48 + (2N - 20)e + (2N + 6)e 2 + 2e 3 , 

a (e,JV) = 2iV(JV + l)(£ + 2iV + 5), ai(e, JV) = (JV + l)(e 2 + 2eiV + 5e + 4N + 12) , (21) 
a 2 (e, JV) = (e - N - 4)(e + 2iV + 3)(e + 2N + 6) 

whic h holds for all JV > 1. T his task can be acc om plished for insta nce by the pack- 
ages Paule and Schorn ( 1995 ). Wegschaider ( 19971) or Schneider ( 2007 ) which are based 
on the cr eative telesc o ping p aradigm presented in IZeilbergerl(jl990ah~ or the paradigm pre- 
sented in iFasenmverl (|1945f ) . Then together with the first two initial values for JV = 1,2, 



S(e,l) 



and 



S(e,2) 



6 



1 



1 



1 



-5 



+ 0(e 3 ), (22) 



e + 6 6 36 

we will be able to compute, e.g., the sum representations of the first 2 coefficients 



F (N) 



3(2JV 2 +4iV+l) 
2N(N+l)(N+2) 



p /at\ 10W < + 52AT- i +63AT+10 
r l\ ly ) — 8AT(N+l)(«+2) 2 



3(-l)~ 
2N(N+l)(N+2) ' 

35i(jV) 



35_i(Af) 



(-l) w (jV-10) 
8N(N+l)(N+2) 2 ' 



2Af(JV+2) ^ 2AT(N+2 

of the e-expansion (14) with t = 0; for more details see Examples 3 and 4 



(23) 
(24) 



In Subsection 3.1 we will develop a recurrence solver which finds the representation 
of the Fi(N) from (14) in terms of indefinite nested sums and products whenever this is 
possible. Afterwards, we combine all these methods to prove Theorem 1 in Subsection 3.2. 

3.1. A recurrence solver for e-expansions 

Restricting the O- notation to formal Laurent series / = X)i*lr f i£l an< ^ 9 ~ Si*l s 9i £l > 
the notation f = g + 0(e t ) for some t E Z means that the order of / — g is larger or 
equal to t, i.e., / — g = Yl^Lt hi £l - Subsequently, K denotes a field with QCKin which 
the usual operations can be computed. We start with the following 

Lemma 1. Let \i € N, and let <Xo(e, JV), . . . , a<j(e, N) G K[e, N] be such that 0^(0, k) ^ 
for all k G N with k > fi. Let h t , ■ . . , h u : N —> K (t, u G Z with t < u) be functions, 
and let Cj ; fc G IK with t < i < u and \i < fc < /i + d. Then there are unique functions 
.Ft, . . . , F u : N — > K (up to the first /i evaluation points) such that Fi(k) = Ci t k for all 
t < i < u and \i < k < d + [i and such that for T(e, JV) = Y^=t -fi(-V) £i we have 

a (e, N)T(e, JV) + • ■ • + a d (e, N)T(e, N + d) = h (N) + ■ ■ ■ + h u (N)e u + 0(e u+1 ) (25) 

for all JV > fi. If the Jii(JV) are computable, the values of the Fi(N) with N > /j, can be 
computed by recurrence relations. 

Proof. Plugging the ansatz T(e, JV) = Y^i=t ^(JV)£ l into (25) and doing coefficient com- 
parison w.r.t. e yields the constraint 



a (0,JV)F t (JV) 



■a d (0,JV)F t (JV + d) = h t {N). 



(26) 



Since a,d(0, JV) is non-zero for any integer evaluation JV > fi, the function Fq : N — > K is 
uniquely determined by the initial values F t (fi) — Ct^, ■ ■ ■ , F t (fi + d — 1) = Ct^+d-i ~ up 
to the first fi evaluation points; in particular the values F t {k) for k > fi can be computed 
by the recurrence relation (26). Moving the F t (N)e t in (25) to the right hand side gives 
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a (e,N) Fi(N)e* + ■ ■ ■ + a d (e, N) ^ F^N + d)e l 

i=t+l »=f+l u 

a (e, N)h t {Ny + ■■■ + a d (s, N)h t (N + d)e*] + ^ 



denote the coefficient of e l on the right side by hi. Since the coefficient of e* on the left 
side is 0, it is also on the right side and we can write 

u u u 

a (e,N) J2 F l {N)e l + --- + a d (e,N) J2 F^N + d)e l = £ h t {N)e l + 0(e u+1 ) 

i=t+l i=t+l i=t+l 

for all TV G N with N > fx. Repeating this process proves the lemma. □ 

Example 3. Consider the recurrence (20) with the coefficients (21). Then by Lemma 1 
there are unique functions F (N) and Fi(JV) with T(N) = F (N) + F e (N)e such that 
T(e, 1) = 2, T(e,2) = 1 + \e and 

ao(e, N)T(e, N) + oi(e, JV)T( £ , JV + 1) + a 2 (e, JV)T( £ , A^ + 2) = fc(e, JV) + 0(e 2 ) (27) 

hold for JV > 1. In particular, by setting e = 0, we get 

a (0,iV) J F 1 o(iV)+ai(0,JV)Fo(JV + l) + a2(0,JV)Fo(JV + 2) = -24JV - 48; (28) 

the values of Fo(JV) can be computed with (28) and the initial values Fo(l) = 2, Fq(2) = 1. 

At th is point we exploit algo rithms from IPetkovsekl (1992); Abramov and Petkovsek 
(jl994l ) ; ISchneiderl (|200ll l2005bl ) which can constructively decide if a solution with certain 
initial values is expressible in terms of indefinite nested products and sums. To be more 
precise, with the algorithms implemented in Sigma one can solve the following problem. 

Problem RS: Recurrence Solver for indefinite nested product-sum expressions. 
Given a (N), a d (N) £ K[N]; given (i£N such that a d (k) / for all k £ N with N > fi; 
given an expression h(N) in terms of indefinite nested product-sum expressions which can be 
evaluated for all JV £ N with JV > /x; given the initial values (c M , . . . , c M _|_d_i) which produce the 
sequence (ci)i>p G K N by the defining recurrence relation 

a (N)c N + ai(JV)cjv+i + ■ ■ ■ + a d (N)c N+d = h(N) VJV > fi. 

Find, if possible, A G N with A > /j, and an indefinite nested product-sum expression g(N) such 
that g(k) = ct for all k > A. 

Remark. Later, we will give further details only for a special case that occurred in almost 
all instances of our computations related to Fcynman integrals; see Theorem 3. 

Example 4. With the input F (l) = 2,F (2) = 1 and (28) Sigma computes the so- 
lution (23). Plugging this partial solution T(e,N) = F (N) + ... into (27) and doing 
coefficient comparison leads to 

V" mp(jv _ -10JV 4 - 98JV 3 - 344JV 2 - 511JV - 267 3{-l) N (3N + 7) 

^^(0,^)^1(^+1) - (7V + 2)(iV + 3)(jV + 4) (AT + 2)(AT + 3)(JV + 4)' 

Then together with Fi(l) = 0, F 1 (2) = 1/6, Sigma finds (24). Since also (2) satisfies (27) 
with the same initial values (22), the first two coefficients of the expansion of (2) are 
equal to Fq(N) and Fi(N) by Lemma 1. 
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This iterative procedure can be summarized as follows. 

Algorithm FLSR (Formal Laurent Series solutions of linear Recurrences) 

Input: /Li G N; oo(e, N), ... , Od(e, iV) G K[s, AT] such that a d (0, k) ^ for all G N with fe > /x; 

indefinite nested product-sum expressions ht(N), . . . , h u (N) [t, u € Z with t < u) which can be 

evaluated for all N G N with N > fj,; dj G K with t < i < u and < j < fx + d. 

Output (r, X,T(N)): The maximal number r G {t — 1, 0, . . . , u} s.t. for the unique solution 

T ( N ) = T,i=t F i( N ) ei with F i( k ) = c i,k for all /i < < ^ + d and with the relation (25) 

the following holds: there are indefinite nested product-sum expressions that are equal to 

Ft(N), . . . , F r (N) for all N > A for some A > fi; if r > 0, return such an expression T(N) 

for T(N) together with A. 

(1) (Preprocessing) By Lemma 1 we can compute as many initial values a t k '■= Fi(k) for 
k > /i as needed for the steps given below (at most A — fi extra values are needed). 

(2) Set r := t, X := fx, and f(N) := 0. 

(3) Note that (F r (iV))]v> M is defined by the initial values F r (N) (A < N < d + A) and the 
recurrence 

ao(Q,N)F r (N) + --- + a d (0,N)F r (N + d) = h r (N) (29) 
for all N G N with N > A; see the proofs of Lemma 1 or Theorem 2. By solving problem RS 
decide constructively if there is a A' > A such that F r (N) can be computed in terms of 
an indefinite nested product-sum expression F r (N) for all N G N with N > X' . 

(4) If this fails, RETURN (r - 1, A, T(N)). Otherwise, set f(N) := f(N) + F r (N)e r . 

(5) If r = u, RETURN (r, A, f(N)). 

(6) Collect the coefficients (product-sum expressions) w.r.t. e l for all i (r + 1 < i < u): 

u 

h'i(N) := coeff (-\a (e, N)F r (N) + ■■■ + a d (e,N)F r (N + d)j + ^ ^(iV)e l ,e l ). 

i — r + l 

(7) Set hi := /i- for all r + 1 < i < u, set r := r + 1 and GOTO Step 3. 

Theorem 2. The algorithm terminates and fulfills the input-output specification. 

Proof. We show that entering the rth iteration of the loop (r > t) we have for all N > X 
that 

u u u 

a (e, N) Fi(Ny + ■■■ + a d {e, N) ^ Fi(N + d)e l = £ h(Ny + 0(e u+1 ) (30) 

i—r i—r i—r 

where the h r (N), . . . , h u (N) are given explicitly in terms of indefinite nested product- 
sum expressions. Moreover, we show that the obtained expression T(N) = Y^iZt Fi(N)e l 
equals the values Xa=t Fi{N)e % for each N > X. For r = t this holds by assumption. Now 
suppose that these properties hold when entering the rth iteration of the loop (r > t). 
Then coefficient comparison in (30) w.r.t. e r yields the constraint (29) for all N > X as 
claimed in Step 3 of the algorithm. Solving problem RS decides constructively if there 
is a A' > such that F r (N) can be computed by an expression in terms of indefinite 
nested product-sum expressions, say F r (N), for all N with N > A'. If this fails, F r (N) 
cannot be represented with such an expression and the output (r — l,A,T(iV)) with 
T(N) = J2i=t Fi(N) is correct. Otherwise, the indefinite nested product-sum expressions 
Fi(N) for t < i <r give the values Fi(N) for all N e N with N > A'. Now move the 
term F r (N)e r in (30) to the right hand side and replace it with F r (N)e r . This gives 

u u d 

a (e,N) F(N)e l + --- + a d (e,N) ^ F(N + d)e l = - ]T 0i ( e , N)F r (N + i) 

i—r+l i—r-\-l i—0 



12 



+ + 0{e u+1 ) =: h r+ i(N)e r+1 + ■■■ + h u {N)e u + 0(e u+1 ) 

i—r 

for all N > A' where h r+ i(N), . . . , h u (N) are given in terms of indefinite nested product- 
sum expressions that can be evaluated for all N £ N with N > A'. By redefining the 
hi(N) as in Step 7 of the algorithm we obtain the relation (30) for the case r + 1. □ 

Algorithm FLSR has been implemented within the summation package Sigma. E.g., the 
expansion for the sum (19) with s = 0, t = 1 and start = 1 is computed by 

GenerateExpansion[a (e,7V)S'[Ar] + ai(e, N)S[N + 1] + a 2 (e,N)S[N + 2], 

{-2AN - 48, 2N - 20}, S[N], {e, s, t}, {start, {{2, 1}, {0, 1/6}}}]; 

here the a^e, N) stand for the polynomials (21), {-24A - 48, 2N - 20} is the list of the 
first coefficients on the right hand side of (20), and start tells the procedure that the 
list of initial values {{2, 1}, {0, 1/6}} from (22) corresponds to N = 1, 2. 

As demonstrated already in Example 4 the following application is immediate. 

Corollary 1. For each nonnegative N, let S(e, N) be an analytic function in e through- 
out an annular region centered at with the Laurent expansion S(e, N) = J2iZt fi(-^) £% 
for some t <G Z, and suppose that S(e, N) satisfies the recurrence (25) with coefficients and 
inhomogeneous part as stated in Algorithm FLSR for some fi £ N; define c^fe := Fi(k) for 
t < i < u and fi < k < fi + d. Let (r, A, Yll=t Fi(N)£ l ) be the output of Algorithm FLSR. 
Then f t (k) = F^k) for all t < i < r and all k £ N with k > A. 

For further considerations we restrict to the following special case. We observed -to our 
surprise- in almost all examples arising from Feynman integrals that the operator 

d 

J2 ai(0, N)Slf = c(N)(S N - b d (N))(S N - b d ^(N)) ...(S N ~ h(N)) (31) 

i=0 

with the shift operator Sn factorizes completely for some b-] .. . . ,b d ,c € K(AQ ; the ra- 



tional functions can be computed by Petkovsek's algorithm Petkovsek ( 1992f) . In this 



particular instance we can construct immediately the complete solution space of 

a (0, N)F(N) + ■■■ + a d (0, N)F(N + d) = X(N) (32) 

for a generic sequence X(N). Namely, choose /ij G N such that the numerator and 
denominator polynomial of bi(j) have no zeros for all evaluations j € N with j > 
and take A := maxi<i<d/i; + 1. Now define for 1 < i < d the hypcrgcomctric terms 
hi(N) = II^La ' — !)■ Then bv lAbramov and Petkovsekl ( 19941 ) one gets the d linearly 
independent solutions 

N-l , ,. > id-2-l , ,. . 

H 1 (N):=h 1 (N),...,H d (N):=h 1 (N)Y,irr^--- E h T I U ^ 

^hx(ix + l) . ^ h d {td-i + 1) 

of the homogeneous version of (32), and the particular solution 

P(N) = /ll< ' iV - > V 1 '"y 1 fedfa-i) '"v^ 1 X(i d ) 

1 '■ c(N) 2^ hl { il + i)"\2^ hd _ l{ld _ 1 + l) 2^ h d (i d + l) 1 1 

»1=A z d _i=A 2d = A 
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of (32) itself. In other words, the solution space of (32) is explicitly given by 

{c 1 H 1 (N) + --- + c d H d (N) + P(N)\c 1 ...,c d eK}; (35) 

here the nesting depth (counting the nested sums) of Hi is i — 1 and of P is d. 
Given this explicit solution space (35) we end up with the following result. 

Theorem 3. Let h t (N), h t+ i(N), . . . with t £ Z be functions that are computable 
in terms of indefinite nested product-sum expressions where the nesting depth of the 
summation quantifiers of hi(N) is di\ let ai(e,N) £ K[e,iV] be such that the operator 
factors as in (31) for some c, hi G ~K(N), c ^ 0. If S(e, N) = YliLt Fi(N)e l is a solution of 

a (e, N)S(e, N) + ■ ■ ■ + a d (e, N)S(e, N + d) = ht(N)e* + h t+1 (N)e t+1 + . . . , (36) 

for some functions Fi (N) , then the values of Fi (TV) can be computed by indefinite nested 
product-sum expressions Fi (N) . The depth of the ^ (iV) is < ma,x t < j<i(dj + (i—j + l)d)). 

Proof. Choose with fi > d such that a d (k) ^ for all integers k > /i and such that 

the sequences hi(k) can be computed for indefinite nested product-sum expressions for 
each k > fi. Consider the rth iteration of the loop of Algorithm FLSR. Since F r (N) is a 
solution of (32) with X(N) = h r (N) for all N > 7, F r (N) is a linear combination of (35). 
Taking the first d initial values F r (/j,), . . . , F r (fj, + d — 1) the Ci are uniquely determined. 
Induction onrgN proves the theorem. The bound on the depth is immediate. □ 

If the operator (29) factorizes as stated in (31), Alg. FLSR can be simplified as follows. 

Simplification 1. The factorization (31) needs to be computed only once and the solutions 
Fi(N) can be obtained in terms of indefinite nested product-sum expressions by simply 
plugging in the results of the previous steps. E.g., for our running example, we get the 



generic solution N N 

E 



(-i) , i(2» 1 +i) (-i)'i(2n+i) (-1) 

' ^ ii(n+l) ^ (2i 2 -l) (2i 2 +l) 

Cl l 1 — 1 t 1 — 1 v / 22 — 1 



7V(AT + 2) 1 C2 2iV(ivT2) 2A^VT2) 1 ' 

of the recurrence a (0, N)F(N) + ai(0, 2V)F(iV + 1) + a 2 (0, AO^W + 2) = X(iV) where 
the coefficients are defined as in (21). In this way, one gets the solution Fq(N) in terms 
of a double sum by setting c\ = c 2 = and X(i 2 ) = —24i 2 + 48 in (37), i.e., 

F (N ) = 1 V (zffil±iil) V '(- 1 )' 224 ^ r 38 i 
2N(N + 2)^> h(l + h) A; (-l + 2i 2 ) (1 + 2*2)" 1 j 

One step further, one gets the solution F\(N) in terms of a quadruple sum by setting 
ci = C2 = and plugging the double sum expression 

X(i 2 ) = 2i 2 - 20 - cocff(a (e, 12)^2) + a x {e, i 2 )F (i 2 + 1) + a 2 (e, i 2 )F (i2 + 2), e) 

into (37). Similarly, one obtains a sum expressions of F 2 (N) with nesting depth 6. 

Minimizing the nesting depth. Given such highly nested sum expressions, the summation 
package Sigma finds alternative sum representations with minimal nesting dept h. The un 



deriving alg orithms are based on a refined difference field theory worked out inlSchneider 
( 2008 , 2010 ) that is adapted from Karr's ITS- fields originally introduced in Karr ( 198ll ). 
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E.g., with this machinery, we simplify the double sum (38) to (23), and we reduce the 
quadruple sum expression for F±(N) to expressions in terms of single sums (24). 

Simplification 2: The solutions (33) of the homogeneous version of the recurrence (32) can 
be pre-simplified to expressions with minimal nesting dept h by the algorithms mentioned 
above. Moreover, using the a l gorith mic theory described in lKauers and Schneiderl (l2006h 



the algorithms in ISchneiderl (|2008l ) can be carried over to the sum expressions like (34) 
involving an unspecified sequence X{id). With this machinery, (37) simplifies to 

Cl C 2 ( — 1J T Zui 1 = l (2i 1 -l)(2i 1 +l) V L ) l^i 1 = l (2i 1 -l)(2i 1 +l) 



N(N + 2) 2N(N + 1)(N + 2) 2N(N + 2) 2N(N + 1) (N + 2) 

Performing this extra simplification, the blow up of the nesting depth for the solutions 
F (N) 7 Fi(N), F 2 (N), . . . reduces considerably: instead of nesting depth 2, 4, 6, . . . we get 
the nesting depths 1, 2, 3, .... In particular, given these representations the simplification 
to expressions with optimal nesting depth in Step 2 also speeds up. 

For simplicity we assumed that the a,i(s,N) are polynomials in e. However, all ar- 
guments can be carried over immediately to the situation where the ai(e,N) are formal 
power series with the first coefficients given explicitly. Moreover, our algorithm is applica- 
ble for more general sequences cii(N) and hi(N) whenever there are algorithms available 
that solve problem RS. E.g., if the coefficients a,i(N) itself are express ible in terms of 
i ndefin ite nested product-sum expression, problem RS can be solved by lAbramov et al 



(|201ll ). and hence Algorithm FLSR is executable. 



3.2. An effective method for multi-sums 

For a multi-sum S(e,N) with the properties (l)-(5) from Assumption 1 and with the 
assumption that it has a series expansion (14) for all N > A for some A £ N, the ideas of 
the previous section can be carried over as follows. 



Step 1: Fi nding a recurrence. B y WZ-theory (jWilf and Zeilbergerl . ll992l Cor. 3.3) and ideas 



given in (jWegschaiderl 119971 Theorem 3.6) it is guaranteed that there is a recurrence of 
the form 

a (e, N)S(e, N) + ■ ■ ■ + a d {e, N)S(e, N + d)=0 (39) 
with coefficients a,i(e, N) £ K[e, N] for the multi-sum S(e, N) in N that can be computed, 
e.g., by Wegschaider's algorithm; for infinite sums similar arguments have to be applied 
as in Step 2.2 of Section 4. Given such a recurrence, let fi £ N with /i > A such that 
a d (Q, N)^0 for all N £ N with N > pi. 

Step 2: Determining initial values. If the sum (13) contains no infinite sums, i.e., s = 0, 
the initial values Fi(k) in <S(e, k) = Y^Lt Fi{^) £l for fc = /i, /.t + 1, . . . can be computed 
immediately and can be expressed usually in terms of rational numbers. However, if 
infinite sums occur, it is not so obvious to which values these infinite sums evaluate 
for our general input class- by assumption we only know that the Fi{k) for a specific 
integer k > fi are real numbers. At this point we emphasize that our approach works 
regardless of whether we express these sums in terms of well known constants or we 
just keep the symbolic form in terms of infinite sums. In a nutshell, if we do not know 
how to represent these values in a better way, we keep the sum representation. However, 
whenever possible it is desirable to rewrite these sums in terms of known values or special 
functions. E xamples are harmonic sums which are known as l imits for the external index 
N — > oo, see Blumlein and Kurth ( 19991 ): Vermaserenl ( 19991 ). to yield Euler-Zagier and 
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multiple zeta value s, cf. iBliimlein et al. (l2010h and references therein, and generalized 



harmonic sums, see Moch et al.l ( 2002 ) which give special values of S'-sums. In massive 2 



loop computations and for the simpler 3-loop topologies these are th e only known classe s, 
whereas extensions are known in case of more massive lines, cf. e.g. Broadhurst ( 19991 ). 

Step 3: Recurrence solving. Given such a recurrence (39) together with the initial values 
of S(e, N) (hopefully in a nice closed form) we can activate Algorithm FLSR. Then by 
Corollary 1, we have a procedure that decides if the first coefficients of the expansion are 
expressible in terms of indefinite nested product-sum expressions. 

Summarizing, we obtain Theorem 1 stated already in the beginning of this section. 
As mentioned already in the introduction, the proposed algorithm (see steps 1,2,3 from 
above) is not feasible for our examples arising form particle physics: forcing Wegschaider's 
implementation to find a homogeneous recurrence is extremely expensive and usually fails 
due to the insufficient computational resources. Subsequently, we relax this restriction 
and search for recurrence relations which are not necessarily homogeneous. 



4. Finding recurrence relations for multi-sums 

Given a multi-sum S(N) of the form (13) we present a general method to compute a 
linear recurrence of S(N). Here the challenge is to deal with infinite sums and summands 
which arc not well defined outside the summation range. We proceed as follows. 

Step 1: Finding a s ummand recurrence. The sum (13) fits the input clas s of the algorithm 



Wegschaiderl (| 1997T) , an extension of multivariate WZ-summation due to lWilf and Zeilberger 
(119921 ). This allows us to compute a recurrence for the hypcrgcometric summand of (13). 
Before giving further details, we recall that an expression T (N, a, j, e) is called hyper- 
geometric in N 7 a 7 j 7 if there are rational functions r VifliV {N,a,j,e) £ K(N,a,j,e) such 
that r(N+i N J+t!j+r,.£) = r v,ii,r,{N,(T,3,e) at the points (v,n,T)) & Z r + S + 2 where this ra - 
tio is defined. Then the Mathematica package MultiSum described in IWegschaiden (|l997t ) 



solves the following problem by coefficient comparison and solving the underlying system 
of linear equations. 

Given a hypergeometric term T (TV, a, j, e), a finite structure set § C N s+r+2 (w.l.o.g. 
we restrict to positive shifts) and degree bounds B G N, /3 <G N s , b £ W +1 . 
Find, if possible, a recurrence of the form 

c u , v , w (N,a,j,e)J 7 {N + u,cr + v,j + w,e) = (40) 

with polynomial coefficients c u ^ ViW S M\N,a,j,e], not all zero, where the degrees of the 
variables N, jt and er; are bounded by B, (3i and hi, respectively. 

Remark 3. (1) In general, choosing S large enough, there always exists a s ummand re- 



curren ce (40) for proper hypcrgcometric summands T (see Footnote 1) due to lWilf and Zeilberger 



(1992). In all our computations we found such a recurrence by setting the degree bounds 
to 1, i.e., B = Pi — hi = 1. 

(2) To determine a small structure set § C N s + r + 2 which provides a solution w.r.t. our 
fixed degree bounds, A. Riese and B. Zimmermann enhanced the package MultiSum by 
a method based on modular computations. In this way one can loop through possible 



16 



choices inexpensively until one succeeds to find such a recurrence (40). 

Next, the algorithm successively divides the polynomial recurrence operator (40) by 
all forward-shift difference operators 

A ai F(N, a, j, e) := T (N, a x , . . . , o l + 1, . . . , a s , j, e) - J"(7V, a, j, e) 

for 1 < i < s, as well as by similar A-operators defined for the variables from j% which 
have finite summation bounds. 

At last we obtain an operator free of shifts in the summation variables (cr, j) called 
the principal part of the recurrence (40) which equals the sum of all delta parts in the 
summation variables from (cr, j), i.e., 



m, cr, j +n,e) 

meS' 1=0 ^ (m,ra)eS; ' 

+ b m ^ n {N,aJ,e)T{N + m,a + k,j + n,e)j (41) 

»=1 ^(m,fc,n)eSi ' 



where the coefficients a m , usually not all zero (see Remark 4.2), b m ^, n and d m ^ n are 
polynomials and the sets S' C N, §i C anc j gj ^ f$ r + 2 are finite. Recurrences of 

the form (41) satisfied by the hypergeometric summand are called certificate recurrences 
and have polynomial coefficients a m (s,N) free of the summation variables from (cr, j), 
while the coefficients of the delta-parts are polynomials involving all variables. 

Remark 4. (1) In principle, the degrees of the polynomials b m} k, n and d m . n arising 
in (41) can be chosen arbitrarily large w.r.t. oi and j'j. However, in Step 2 we will 
sum (41) over the input range and hence we have to guarantee that the resulting sums 
over (41) are well defined. As a consequence, the degrees of the d m , n and b„ lt k,n w.r.t. 
the variables <Ji have to be chosen carefully if in (15) one of the constants Ci is zero. 
As mentioned earlier, for such situations we restrict ourselves to the case s = 1. In this 
case, the degree in the b mi k, n should be smaller than the constant d\ from (15) and the 
degree in the d m , n should be not bigger than the constant o from ( 17). To control this 
total bound b := min(<ii — 1, o), we exploit the following observation ( Wegschaiderl fl997l 



p. 43): While transforming (40) to (41) by dividing through the operators (4), one only 
has to perform a simple sequence of additions of the occurring coefficients in (40), and 
thus the degrees w.r.t. the variables do not increase. Summarizing, if we choose (3\ in our 
ansatz such that /3i < &, the degrees in the b m ^^ n and d m ^ n w.r.t. the variable o\ arc 
smaller than b. 

(2) In general, it might happen that t he principal part is 0, i.e., we get a trivial remain- 



der within the operator divisions. In (|Wegschaiderl . 119971 Thm. 3.2) this situation was 



resolved at the cost of increasing the degrees w.r.t. some of the variables. If within this 
construction the degree w.r.t. cri increases too much, manual adjustment is needed (e.g., 
force the structure set to be different or change the degree bounds manually). However, 
this exotic case never occurred within our computations. 
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Example 5. For the sum 



V Jl + 1 I ( 4 - £ )A+ji(§ +4) Jo+Jl 



with the discrete parameter AT > 3 and e > the package MultiSum computes the 
summand recurrence 

(e - 2N)NF(N,j ,j 1 ) - (e -N - 3)(e + 2iV + 2)^(7V + 1, j ,Ji) 

= A J0 [(e 2 + j e + e - 2 3l - 2j N - ifrN - 12N - 6)F(N + 1, j , n )} 

+ A h [(e - 2N)V + h -N + 1)T(N, jo, h) 
+ {-2N 2 + eN + 2j N + A n N + AN-2e- ej + 2j l )F{N + 1, j , ji]). (43) 

Step 2: A recurrence for the sum. Taking as input the certificate recurrences (41) we 
algorithmically find the inhomogcneous part of the recurrence satisfied by the sum (13) 
which will contain special instances of the original multi-sum of lower nesting depth. 

The recurrence for the multi-sum (13) is obtained by summing the certificate recur- 
rence (41) over all variables from (<r, j) in the given summation range 1Z C Z s+r+1 . Since 
it can be easily checked whether the summand J- satisfies the (41), the certificate re- 
currence also provides an algorithmic proof of the recurrence for the multi-sum S(N,e). 
In particular, since we set up the degrees of the coefficients in (41) w.r.t. the variables 
accordingly, see Remark 4, it follows that the resulting sums arc analytically well defined. 

To pass from the certificate recurrence to a homogeneous or inhomogcneous recur- 
rences for the sum, special emphasis has to be put on the A-operators. In particular, the 
finite summation bounds appearing in (13) lead to an inhomogcneous right hand side af- 
ter summing over the summand recurrence (41). A method to se t up the inh omogeneous 
recurrences for the summation problems (13) was introduced in (|Stanl . l2010l Chapter 3). 



We summarize the steps of this approach implemented in the package FSums. 
In this context, we use tuples to denote multi-dimensional intervals. The range repre- 
sented by the tuple interval [i,k] is the Cartesian product of the intervals defined by 
the components i, k € Z n . More precisely, [i,k] := [ii,/ci] x [12,^2] x ■•• x [i n , k n ] where 
[ij, kj] = {ij, + . . . , kj}. Often when working with nested sums, summation ranges for 
inner sums will depend on the value of a variable for an outer sum. Intervals whose end- 
points are defined by tuples are not enough to represent the summation ranges for these 
sums. We will use a variant of the cartesian product notation to denote such a summation 
range. Namely, to refer to a variable associated to a range, we will specify it as a subscript 
at the corresponding interval and use x signs instead of the x symbols. For example, the 
range for the sum (18) can be written as [0, 00) x [0, AT — 3] J0 X [0, A^ — jo — 3] k [0, jo + 1] ■ 
We also introduce this notation for the initial range of the sum (13) as 

K — KaXfcj (44) 

where lZ a := \p, 00) and IZj = [qo, L\(N)] X • • • x [q r , L r (N,jo, . . . , jV-i)], are the infinite 
and the finite range, respectively. 

Step 2.1: Refining the input sum. As indicated earlier, we consider the summands from (13) 
as well-defined only inside the initial input range TZ C Djr where Djr denotes the set of 
well-defined values for the proper hypcrgcomctric function T . Because of this restriction 
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we need to determine a possible smaller summation range over which we are allowed to 
sum the certificate recurrences (41). 

Example 6. We illustrate this phenomenon by our concrete example (42). Let us start 
by summing over the initial summation range 1Z = [0, N — 3] J0 k [0,N — 3 — jo] over 
the delta parts on the right hand side of the recurrence (43) which is of the form (41). 
For this we denote the polynomial coefficients inside the delta parts A J0 and Aj 1 with 
e(N, jo, ji, e) and d\{N,jo,ji,e), d.2(N,jo,ji,e), respectively. By summing over the first 
term inside the A^-part and using the telescoping property, we have 



N-3 N-3-j N-3 

E v 



^ Aj^iNJoJ^e^iNJoj!)] = ^ (d x (iV,io,ii,£)J-(iV,io,ii)) 



ji=N-2-j 
31=0 



30=0 %°3 J °=° N-3 

= J2di(N,j ,N-2- jo , e)J 7 (N, j , N - 2 - j ) - ^ di (N, j ,0, e)F(N, j , 0) 

3o=0 jo=0 

where we use the short-hand notation J2k=o -F(k, : = J2k=o -^"(^i B)—J2k=o -^C ' A)- 

We observe that, after telescoping, the upper bound N — 2 — jo for ji translates into 
a term outside the original summation range. To work under the assumption that our 
summand T(N, jo, ji) is well-defined only inside its range 72., we need to adjust the range 
over which we sum the certific ate recurren ce or shift this relation with respect to the free 
parameter N. As discussed in (Stan, 20101 Chapter 3), the approach based on computing 



a smaller admissible summation range is more efficient since it leads to fewer new sums 
in the inhomogencous parts of the recurrences. 

In the case of our example S(e, N), we consider the new range TV = [0, N — 4] J0 x [ 0, N 



7n — 4 1 . As a consequence we compute separately a single sum which was called in ([Stan 



20101 Chapter 3) a sore spot, 



N-iN-4-jo N-3 

S(e,N) = J2 J2 W io.Ji) + HN,jo,N-j Q - 3). (45) 

30=0 3i=0 3o=0 

In general, the package FSums contains an algorithm that determines the inevitable 
summation range and computes the necessary sore spots for sums of the form (13); 
these extra sums with lower nesting depth have to be considered separately (see also the 
DIVIDE step in our method described in Section 5). Subsequently, we denote the sum 
over the restricted range 1Z' by S'(e, N). 

Step 2.2: Determining the inhomogeneous part of the recurrence. Summing a certificate 
recurrence of the form (41) over the restricted range 1Z' determined in the previous 
step leads to a recurrence for the new sum S'(e,N). The inhomogeneous part contains 
special instances of this sum of lower nesting depth. Next, we introduce the types of sums 
appearing on the right hand side. 

Step 2.2.1: The finite summation bounds. Shift compensating sums are the first side-effect 
of nonstandard summation bounds. They appear when we sum over the left hand side of 
the recurrence over a given definite range, because our upper summation bounds depend 
on the other summation parameters. 

Example 7. Subsequently, we will illustrate these aspects with our running example 



19 



(42). As deduced from Step 2.1, we continue from now on with the new sum 

7V-4 N-4-jo 

S'(e,N)=J2 E Wio.Ji). (46) 

jo=0 3 i=0 

When we sum the certificate recurrence (43) over the restricted range 1Z\ we obtain 

JV-4JV-4-J0 N-3 

J2 E + 1 'J'o,Ji) =S'(e,N + 1) -J2 HN+1,3, N-3- j). (47) 

jo=0 ji=0 j=0 

Compensating sums of this form appear only in the case of upper summation bounds 
depending on the free variable N . After summing over the left hand side of the re- 
currence, we will move the resulting compensating sums, with a change of sign, to the 
inhomogeneous part. 

Example 8. Including the new shifted sum as the first term of the output, the following 
procedure of FSum delivers the right hand side of (47) 

ln[l]:= ShiftCompensatingSums[F[AT,jo,ii],{{io,0, N - 4}, {ji , 0, N - 4 - jo}}, AT, 1] 

Out[i]= SUM[7V + 1] + FSum[-F[l + N,jo, -3 - jo + N], {{j ,0, -3 + TV}}]. 

Note that we use the structure FSum to store sums with nonstandard boundary conditions 
of the form (13). This data type contains two components, the summand and a list 
structure for the summation range. The nested range is stored in the order given in (13), 
starting with the infinite sums and ending with the sums with finite summation bounds 
in the order of their dependence. 

When summing over the A-parts we generate two types of sums on the right side of 
the recurrence, the A-boundary sums and the so-called telescoping compensating sums. 

Example 9. When summing over the A JO -part of the recurrence (43), we get 

JV-3 N-3-ja 

E ^io[<N,j ,ji,e)r(N + l,jo,ji)] 
io=o 3i=0 

JV-2 JV-2-io N-3 N-3-j 

= E E e{N,j ,j 1 ,e)F(N + l,j ,j 1 )-^2 E <N,j ,j 1 ,e)F(N+l,j ,j 1 ). 

30 =1 31=0 jo=0 31=0 

Now one sees that exactly the sum with the summation index jo cancels and one obtains 



JV-3-jo 



J0=N-2 



+ ^e(N,j ,N-2-j ,e)F(N+l,j ,N-2-j ). 



30=0 



30 = 1 



(e(JV ) io,ii,e)W + l,j'o,Ji)) 

31=0 

Because of the structure of the summation bounds for the nested sums (13) we can use 
again our procedure Shif tCompensatingSums to generate the shift compensating sums 
and to read off the telescoping compensating sums. This connection becomes clearer 
when we consider the more involved sum (18) (with its restricted range N — 4 instead of 
its original range N — 3) and apply, e.g., the A J0 -operator: 

i'o=JV-3 

oo JV-4JV-jo-4 jo oo N-jo-i jo 

E E E E a a \hn, ,30,31,32)] = E E E JF ( N > > ^ > * > & ) 

cto=0 j a =0 ji=0 j 2 =0 cr o =0 ji=0 j 2 =0 . 

3o=0 
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oo iV-3j'o-l oo JV— 3JV-j'o-4 

+ F(N,<T ,jo,N-j -3,j 2 ) £ F{N,<roJo,h,joY, 

co=0 jo = l ]2=0 o- =0j'o = l ji=0 

note that the first element on the right side of this identity produces the A-boundary 
sums while the last two are due to telescoping compensation. More precisely, with 

ln[2]:= ShiftCompensatingSums[F[JV, er , jo — 1, ji , jz\ , {{cro, 0, oo}, {j±, 0, N — jo — 4}, 
{j2,0, jo}}/ -jo -»■ (jo - 1), jo, 1] 

Out[2]= {FSnm[F[N,a ,j a ,j 1 ,j 2 ],{Wo,0,OD},{j 1 ,0,N-A-j },{j 2 ,0,jo}}],FSum[F[N,ao,jo,N-3- 
3o,j2],{{ao,0,oo},{j 2 ,0,jo - 1}}], FSum[-F[/V, a ,j ,ji,jo], {Wo,0, oo}, {ji, 0, TV - 4 - j }}]} 

we obtain exactly this result: the delta boundary sums are obtained by evaluating the 
first entry of the output for jo = and jo = N — 3 and the compensating sums result 
by adding the shifted sum [1, TV — 3]j to the range of the o ther terms in the output. A 
detailed description of these computations can be found in (IStanl . l20ld Alg. 4). 



Step 2.2.2: The infinite summation bounds. To sum over the delta parts in (41) coming 
from the summation variables Ui, e.g., A ai b m ^k,n{N, a, j, e)J-(N + m, a + k, j + n, e) we 
have to ensure that lim (Ti _ i . 00 b m ,k,n{N, a, j, e) .F(iV + m,a + k,j + n, e) exists. Looking at 
the asymptotic conditions (15) of the input sum (13), there will be no problem if Cj > 0. 
However, if the constant q is zero, we need to verify that the degrees of the polynomial 
coefficients b m ,k,n appearing in the respective A CTi -part are smaller than the bound Pi. 
As worked out in Remark 4 this property is guaranteed by our ansatz. 

The above sections introduced the types of sums, i.e., shift and telescoping compen- 
sating sums as well as delta boundary sums, which will appear on the right hand side of 
the inhomogeneous recurrences satisfied by summation problems of the form (13) after 
summing over corresponding certificate recurrences (41). A procedure to generate these 
inhomogeneous recurrences is implemented in the package FSums. E.g., the recurrence 
satisfied by the sum S'(e,N), which we denote by SUM [TV], is returned by 

ln[3]:= finalRecS = InhomogenRec [certRecS , {{jo, 0, —4 + N}, {ji, 0, —4 — jo + N}}, N] 

Out[3]= (e - 2JV)JVSUM[JV] + (3 - e + N)(2 + e + 27V)SUM[1 + N] == 

FSum[(l + jo - N)(-e + 2N)F[N,j , 0], {{jo,0, -4 + N}}] + 
FSum[-2(e - 2N)F[N,j , -3 - jo + N],{{j ,0, -4 + N}}] + 
FSum[(e - 2N)(2 + j - N)F[1 + N,jo, 0], {{jo, 0, -4 + N}}] + 
FSum[(6 -e-e 2 + 2ji + 127V + 4ji N)F[l + N, 0, ji], {{j lt 0, -4 + N}}] + 
FSum[(3 - e + N)(2 + e + 2N)F[1 + N, j , -3 - jo + N], {{jo, 0, -3 + N}}] + 
FSum[(e + e 2 + 2j + ejo - 2N + 2j N - 4/V 2 )F[l + N, jo, -3 - jo + N], {{jo, 1, -3 + N}}} + 
FSum[-((6 + 2e + 2jo + ejo + 6/V-e/V + 2jo/V-2/V 2 )F[l + /V,j ,-3-jo + /V]),{{jo,0, -4 + 7V}}]; 

here certRecS stands for the certificate recurrence (43). 



5. An efficient approach to find e-expansions for multi-sums 

Let S(e, N) be a multi-sum of the form (13) with the properties (l)-(5) from Assump- 
tion 1 and assume that <S(e, N) has a series expansion (14) for all N > A for some A 6 N. 
Combining the methods of the previous sections we obtain the following general method 
to compute the first coefficients, say F t (N), . . . , F U (N) of (14). 
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Divide and conquer strategy 

(1) BASE CASE: If S(e, TV) has no summation quantifiers, compute the expansion by 
formulas such as (8) and (9). 

(2) DIVIDE: As worked out in Section 4, compute a recurrence relation 

a {e, N)S{e, N) + ■ ■ ■ + a d {e, N)S(e, TV + d) = h(e, N) (48) 

with polynomial coefficients a,i(e,N) G K[e, TV], a m (e,N) ^ and the right side 
/i(e, N) containing a linear combination of hypcrgcomctric multi-sums each with 
less than s + r + 1 summation quantifiers. Note: In some cases, the sum has to be 
refined and some "sore spots" (again with fewer summation quantifiers) have to be 
treated separately by calling our method again; see Step 2.1 in Section 4. 

(3) CONQUER: Apply the strategy recursively to the simpler sums in h(e,N). This 
results in an expansion of the form 

h(e,N) = ht(N)i* + h t+1 (N)e t+1 + ■■■ + h u {N)e u + 0{e u+1 ); (49) 

if the method fails to find the h t (N) , . . . , h u (N) in terms of indefinite nested 
product-sum expressions, STOP. 

(4) COMBINE: Given (48) witiGD (49), compute, if possible, the F t (N), . . . , F U (N) 
of (14) in terms of nested product-sum expressions by executing Algorithm FLSR. 

We illustrate our method with the double sum (42); internally we transform all the 
objects in terms of T(a;)-functions in order to apply expansion formulas such as (8) 
and (9). First, we compute the summand recurrence given in (43). While computing a 
recurrence for the sum itself, it turns out that we have to refine the summation range, i.e., 
our computation splits into two problems as given in (45). We continue with the refined 
double sum (46) and obtain the inhomogeneous recurrence finalRecS given in 0ut[3]. 
Now we apply recursively our method and compute successively expansions for each of 
the single sums on the right hand side; see also Example 2. Adding all the expansions 
termwise gives the recurrence 

(e - 27V)TV<S'(e, TV) - (e - TV - 3)(e + 27V + 2)<S'(e, TV + 1) = 

18(2TV 6 - 37V 5 - 87V 4 + 137V 3 - 47V + 8) 36(2TV 4 + N 3 - 97V 2 - 27V + 4)(-l) Ar 



+ e 
+ 



(TV - 2)(7V - 1)7V(7V + 1)(7V + 2) (TV — 2) (TV - 1)TV(TV + 1)(TV + 2) 

3(7V 8 - 67V 7 - 327V 6 + 207V 5 + 1517V 4 + 147V 3 - 2007V 2 - 287V + 56) 
(TV - 2) (TV - 1)7V(7V + 1) 2 (7V + 2) 2 
6(2TV 6 + TV 5 — 14TV 4 + 97V 3 + 407V 2 - 227V - 28)(-l) JV 365i(7V) 



2 

+ e 



(7V-2)(7V-l)7V(7V + l) 2 (7V + 2) 2 TV + 1 

gS^TV) 2 6(TV - 5)Si(TV) TV 6 (5TV 3 + 487V 2 + 246TV + 568) 



TV + 1 (TV+1) 2 4(7V~ l)(7V-2)(7V + l) 3 (7V + 2) 3 

9 (TV 4 - TV 3 - 4 TV 2 + 4TV + 8) 18 (2TV 4 + TV 3 — 97V 2 - 27V + 4) (-1)' 



(TV - 2) (TV - 1)7V(7V + 1)(7V + 2) (TV - 2) (TV - 1)TV(7V + 1)(7V + 2) 
3637V 6 + 37207V 5 + 36727V 4 - 52807V 3 - 107127V 2 - 45927V - 128 



S 2 (N) 



+ +- 



47V(TV - 1)(TV - 2) (TV + 1) 3 (TV + 2) 3 



+ 0(s 3 ) 



Together with its first initial value S'(e, 4) = — j^£ — y^e 2 Algorithm FLSR computes 
the series expansion of S'(e,N). Finally, we compute the expansion of the extra sum 



3 Cf. Step 2 of Section 3.2 to see how we deal with the initial values. 
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Ylf =a-F(N,jo, N — 3 — jo) with our method, and adding this result to our previous 
computation leads to the final result 



Ar , 81(iV 2 -3JV + 2) 



3(JV 4 - 13JV 3 - 28JV 2 - 32JV + 24) + 9(JV + 3)Si(N) 



8V 3 (JV + 2) JV(JV + 1){N + 2)_ 

9(/V + 3)Si(/V) 2 3(5JV 3 + 36JV 2 +37JV-18)S'i(JV) 9(JV 2 + 3N + 4)5* 2 (JV) 



+ 



4JV(JV + 1){N + 2) 4JV(JV + l) 2 (JV + 2) 2 4JV 2 (V + 1)(JV + 2) 



5JV b + 17JV 5 + 162JV 4 + 208JV 3 + 592JV 2 + 240JV - 288 



+ 0(0 



32JV 4 (V + 2) 2 

Similarly, we compute, e.g., the first two coefficients of the expansion of the sum (18): 
Ijum 3(-l) iv (jV 2 + 2jV-l)g 1 (iV) 

w(£liV) - jv(jv + i) 9( ~ 1} + jv + 

3(-l) JV (4iV + 3) 3(-l) JV (3iV + 2) S_i(JV) 9, 3{-l) N Si{N) 2 



:[C(2) 



2JV JV 2N> 2 (N + 1) 

(-l) Ar (2AT 4 + 34AT 3 + lOliV 2 + 89iV + 2)5i(JV) 3(-l) JV (4iV 2 + 14iV + 13) 



+ 



2N(N + l) 2 (JV + 2) (7V+l)(7V + 2) 

(-1)^ (-30JV 2 -38JV + l)g 2 (/V) 9(-l) JV (2JV + l)g 3 (iV) , 9 S_ 2 (iV) w 

+ 2JV(JV + 1) JV + JV +9(_1) S2 ' l{N) 

6(-l) JV (3JV + 2)S- 2 (JV)5-i(JV) 6(-l) JV (3JV + 2) S-2,-i(JV)l . 2 ^ 



JV JV 
where C(2) = E£x ^ = ^/6- 



Remark 5. In the following we give further comments on our proposed method and 
provide strategics for using it in the context of the evaluation of Feynman integrals. 

1. A heuristic. The conquer step turns our procedure into a method and not into an 
algorithm. Knowing that there is an expansion of S(e, JV) in terms of indefinite nested 
sums and products and plugging this solution into the left hand side of (48) shows that 
also the right hand side of (48) can be written in terms of indefinite nested product-sum 
expressions. But in our method the right hand side is split into various sub-sums and it 
is not guaranteed that each sum on its own is expressible in terms of indefinite nested 
product-sum expressions - only the combination has this particular form. However, for 
our input class arising from Feynman-integrals this method always worked. 

2. A hybrid version for speed-ups. As it turned out, the bottleneck in our computations 
is the task to compute a recurrence of the form (48) with the MultiSum-package. To 
be more precise, in several cases we succeeded in finding a structure set S with the 
corresponding degree bounds for the polynomial coefficients, but we failed to deter- 
mine the summand recurrence (40) explicitly, since the underlying linear system was 
too large to solve. For such situations, we dropped, e.g., the outermost summation quan- 
tifier, say X)o-i=pi an d searched for a recurrence in ay, in particular the variable JV 
was put in the base field K. In this simpler form, we succeeded in finding a recurrence. 
Next, we computed the initial values (in terms of JV) by using another round of our 
method. With this input, Algorithm FLSR found an expansion with coefficients in terms 
of Ft(N, a), Ft+i(N, a), . . . , F U (N, a). To this end, we applied the infinite sum 



J2 F ^ N ^ ( 5 °) 



to the coeffi cients Fj (JV, a) and si mplified these expressions further by the techniques 



described in lAblinger et al.l (|2011bl ). In various situations, it turned out that this hybrid 
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technique was preferable to computing a pure recurrence in 7V or ju st simplifying the 
expressions (12) by using the methods given in lAblinger et al.l (|2011bl ). 

3. Asymptotic expansions for infinite expressions. As mentioned in Remark 2 wc obtained 
also sums of the form (13) which could be defined only by considering a truncated version 
of the infinite sums. For such cases we computed the coefficients Fi(cr,N) as above and 
considered -instead of (50)- the expressions 'YT a=0 Fi(a 1 N) for large values a. To be 
more precise, we computed asymptotic expansions for all these sums and combined them 
to one asymptotic expansion in a. In this final form all the expressions canceled which 
were not defined when performing a — > oo and we ended up with the correct Fi(N). 

4. Dealing with several infinite sums. In all our computations only a single infinite sum 
arose. In principle, our method works also in the case when there are several such sums. 
However, in order to set up the recurrence in Section 4, we need additional properties 
such as (17) for the multivariate case. If such properties are not available, we propose two 
strategics: 4.1 Drop some (or all) of the infinite sums and proceed as explained in point 
2 of our remark. 4.2 Set up the recurrence with formal sums and expand the sums on the 
right hand side: here one can either use the strategies as described in Step 4 of Section 2 
(in particular, if asymptotic expansions have to be computed), or one can proceed with 
the method of this section whenever the sum is analytically well defined. 



6. Conclusion 



We presented a general framework that enables one to compute the first coefficients 
Fi(N) of the Laurent expansion of a given Feynman parameter integral, whenever the 
Fi(N) are expressible in terms of indefinite nested product-sum expressions. Namely, 
starting from such integrals, we described a symbolic approach to obtain a multi-sum 
representation over hypcrgcomctric terms. Given this representation, we developed sym- 
bolic summation tools to extract these coefficients from its sum representation. In order 
to tackle this problem, Wegschaider's MultiSum package has been enhanced with Stan's 
package FSum that handles sums which do not satisfy finite support conditions. More- 
over, given a recurrence relation of the form (36) together with initial values, we used 
Schneider's recurrence solver that decides constructively, if the first coefficients of the 
formal Laurent series solution are expressible in terms of indefinite nested product-sum 
expressions. 

In order to fit the input class of hypcrgcomctric multi-sum packages, we split the 
sums at the price of possible divergencies. We ov ercame this situation b y combining our 
new methods with other tools described, e.g., in lAblinger et al. (l2011bh : see Remark 5. 
Further analysis of the introduced method should lead to a uniform approach that can 
handle in one stroke also solutions in terms of asymptotic expansions. 

The described summation tools assisted in the task to compute two- and simple 
three-loop diagrams, which occurred i n the calculation of th e massive Wilson coeff i- 



cients for deep-inelastic scattering; s ee lAblinger et al. (|2011dI) : iBliimlein et all (<2006h 



Bierenbaum et al.l ([2007, 2009a, 2008) . We are curious to see whether these new summa- 



tion technologies find their application also in other fields of research. 
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